US COVID Data

Data Prep

us <- read.csv("https://raw.githubusercontent.com/JaclynCoate/6373_Time_Series/master/TermProject/Data/USdaily.7.26.20.csv", header = T, strip.white = T)
us <- transform(us, date = as.Date(as.character(date), "%Y%m%d"))
us <- subset(us, select = -c(states, dateChecked, hospitalized, lastModified, total, posNeg, totalTestResultsIncrease, hash))
us[is.na(us)] <- 0
us = us[order(as.Date(us$date, format = "%Y%m%d")),]
head(us)
##           date positive negative pending hospitalizedCurrently
## 187 2020-01-22        2        0       0                     0
## 186 2020-01-23        2        0       0                     0
## 185 2020-01-24        2        0       0                     0
## 184 2020-01-25        2        0       0                     0
## 183 2020-01-26        2        0       0                     0
## 182 2020-01-27        2        0       0                     0
##     hospitalizedCumulative inIcuCurrently inIcuCumulative
## 187                      0              0               0
## 186                      0              0               0
## 185                      0              0               0
## 184                      0              0               0
## 183                      0              0               0
## 182                      0              0               0
##     onVentilatorCurrently onVentilatorCumulative recovered death
## 187                     0                      0         0     0
## 186                     0                      0         0     0
## 185                     0                      0         0     0
## 184                     0                      0         0     0
## 183                     0                      0         0     0
## 182                     0                      0         0     0
##     totalTestResults deathIncrease hospitalizedIncrease negativeIncrease
## 187                2             0                    0                0
## 186                2             0                    0                0
## 185                2             0                    0                0
## 184                2             0                    0                0
## 183                2             0                    0                0
## 182                2             0                    0                0
##     positiveIncrease
## 187                0
## 186                0
## 185                0
## 184                0
## 183                0
## 182                0

Current Hospitalization Realization

Traits:

  • Heavy wandering behavior
  • What appears to be some noise that could be pseudo-cyclic behavior hidden by the large numbers.
ggplot(data = us, aes(x=date, y=hospitalizedCurrently))+
  geom_line(color="orange")+
  labs(title = "Current COVID Hospitalized Cases US", y = "Thousands", x = "") +
    theme_fivethirtyeight()

Stationarity

It is difficult to assume stationarity for this data due to multiple factors. We are working under the assumption that COVID is a novel virus and cases as well as hospitalizations will eventually return to zero. This being said our current modeling techniques do things such as return to the median or mimic the previously seen trends. Also, we see a heavy wandering trend in both new cases and hospitalization would be dependent on this as well as time. We will review the data and see what, if any, non-stationary components reveal themselves and model the data accordingly.

Forecast Current Hospitalizations

It is difficult to assume stationarity for this data due to multiple factors. We are working under the assumption that COVID is a novel virus and cases as well as hospitalizations will eventually return to zero. This being said our current modeling techniques do things such as return to the median or mimic the previously seen trends. Also, we see a heavy wandering trend in both new cases and hospitalization would be dependent on this as well as time. We will review the data and see what, if any, non-stationary components reveal themselves and model the data accordingly.

Variable Review

When reviewing the variables and the scatter plot from our previous EDA we can see that there are some correlations that were expected, for example currentHospitalizations is correlated with the variables that refelct ICU and Ventilaor patients. These metrics wouldn’t exist without hospitalizations. These variables also cannot help up predict hospitalizations because they after occurances. So we will actually be leveraing variables such as positive and positive increase in order to see if there is some sort of correlation between hospitalized patients and the number of positive cases.

  • Positive Tests Trend
ggplot(data = us, aes(x=date, y=positive))+
  geom_line(color="orange")+
  labs(title = "Total Cases US", y = "Millions", x = "") +
    theme_fivethirtyeight()

  • Positive Increase Trend
ggplot(data = us, aes(x=date, y=positiveIncrease))+
  geom_line(color="orange")+
  labs(title = "Increase in COVID Cases US", y = "Thousands", x = "") +
    theme_fivethirtyeight()

MLR w/ Correlated Errors for Currently Hospitalized Patients

Forecast Independent Variables

  1. Forecast Increase Positive Cases
    1. Evaluation: Slowly dampening ACF and heavy wandering consistent with a (1-B) factor. As well as frequency peaks at 0 and .14. Consistent with (1-B) and seasonailty component of 7.
#a
plotts.sample.wge(us$positiveIncrease)

## $autplt
##  [1] 1.0000000 0.9740745 0.9439116 0.9109437 0.8869924 0.8711719 0.8605261
##  [8] 0.8462935 0.8144161 0.7772852 0.7356766 0.7045767 0.6829609 0.6646632
## [15] 0.6432819 0.6088202 0.5688872 0.5293417 0.5001353 0.4761589 0.4592724
## [22] 0.4426271 0.4148783 0.3789333 0.3458221 0.3182895
## 
## $freq
##  [1] 0.005347594 0.010695187 0.016042781 0.021390374 0.026737968
##  [6] 0.032085561 0.037433155 0.042780749 0.048128342 0.053475936
## [11] 0.058823529 0.064171123 0.069518717 0.074866310 0.080213904
## [16] 0.085561497 0.090909091 0.096256684 0.101604278 0.106951872
## [21] 0.112299465 0.117647059 0.122994652 0.128342246 0.133689840
## [26] 0.139037433 0.144385027 0.149732620 0.155080214 0.160427807
## [31] 0.165775401 0.171122995 0.176470588 0.181818182 0.187165775
## [36] 0.192513369 0.197860963 0.203208556 0.208556150 0.213903743
## [41] 0.219251337 0.224598930 0.229946524 0.235294118 0.240641711
## [46] 0.245989305 0.251336898 0.256684492 0.262032086 0.267379679
## [51] 0.272727273 0.278074866 0.283422460 0.288770053 0.294117647
## [56] 0.299465241 0.304812834 0.310160428 0.315508021 0.320855615
## [61] 0.326203209 0.331550802 0.336898396 0.342245989 0.347593583
## [66] 0.352941176 0.358288770 0.363636364 0.368983957 0.374331551
## [71] 0.379679144 0.385026738 0.390374332 0.395721925 0.401069519
## [76] 0.406417112 0.411764706 0.417112299 0.422459893 0.427807487
## [81] 0.433155080 0.438502674 0.443850267 0.449197861 0.454545455
## [86] 0.459893048 0.465240642 0.470588235 0.475935829 0.481283422
## [91] 0.486631016 0.491978610 0.497326203
## 
## $db
##  [1]  14.3997611  15.8955909   7.7642205   8.7050186   3.8198773
##  [6]  -1.1863576   2.6815084  -0.4667544  -1.0582274  -3.9729372
## [11]  -4.4388886  -5.3739230  -7.0574666  -3.9001118  -5.3732547
## [16]  -5.9203161  -6.0265835  -6.3073466  -8.7427682  -8.5315257
## [21]  -8.1694048  -8.8726978  -6.2228751  -5.6255282  -5.7869081
## [26]  -2.3190255  -0.6895940 -18.2510780 -12.7252753 -17.9895004
## [31] -13.1495229 -17.1809062 -12.9879824 -18.5885318 -14.0950716
## [36] -19.2823521 -16.1802599 -22.2193540 -17.1778910 -15.9864725
## [41] -13.8720977 -14.3606703 -12.3482825 -19.0159734 -12.8299009
## [46] -16.1634117 -19.2402721 -16.9940391 -19.9193164 -24.4462952
## [51] -12.2184745 -14.1347405 -11.1407485 -19.5608640 -22.3406178
## [56] -26.2948840 -15.9219942 -19.7896219 -16.3429699 -16.6545137
## [61] -13.2991307 -20.9465604 -12.5392646 -23.3015328 -21.5282804
## [66] -17.4748203 -26.6390880 -24.5923917 -17.5647166 -19.7282494
## [71] -17.6845420 -23.0825883 -19.0743556 -19.4958696 -19.5138333
## [76] -20.4270518 -18.7343253 -16.8982195 -14.1306252 -14.2054286
## [81] -15.0466299 -13.4474119 -14.7662082 -15.2848488 -14.7767269
## [86] -17.2543647 -17.5524024 -16.7680943 -23.6931393 -19.6888637
## [91] -22.9328710 -17.1805763 -15.7422487
## 
## $dbz
##  [1]  12.1744244  11.8068871  11.1912920  10.3234897   9.1987290
##  [6]   7.8134517   6.1690554   4.2794940   2.1855400  -0.0227779
## [11]  -2.1856123  -4.0930116  -5.5814407  -6.6274820  -7.3078903
## [16]  -7.6988921  -7.8526856  -7.8271468  -7.6919809  -7.5061868
## [21]  -7.3030888  -7.0977786  -6.9055053  -6.7545171  -6.6859448
## [26]  -6.7448940  -6.9710578  -7.3933729  -8.0283548  -8.8797888
## [31]  -9.9376266 -11.1745795 -12.5395304 -13.9486256 -15.2804168
## [36] -16.3913300 -17.1665125 -17.5812779 -17.7097183 -17.6695854
## [41] -17.5637030 -17.4562894 -17.3753921 -17.3239844 -17.2913960
## [46] -17.2628022 -17.2262162 -17.1767997 -17.1186074 -17.0640034
## [51] -17.0309217 -17.0382729 -17.1002705 -17.2209453 -17.3905128
## [56] -17.5857017 -17.7761429 -17.9367144 -18.0606665 -18.1649087
## [61] -18.2832201 -18.4521578 -18.6979308 -19.0286543 -19.4319317
## [66] -19.8765383 -20.3176913 -20.7053075 -20.9929698 -21.1442791
## [71] -21.1356798 -20.9585908 -20.6228760 -20.1586258 -19.6118560
## [76] -19.0349178 -18.4769322 -17.9782467 -17.5690679 -17.2703544
## [81] -17.0951532 -17.0493949 -17.1318377 -17.3332302 -17.6350503
## [86] -18.0085946 -18.4157858 -18.8133902 -19.1613062 -19.4324915
## [91] -19.6189508 -19.7292633 -19.7788462
b. Differencing Data: much more stationary data and have surfaced a seasonaly component = 7 seen on ACF peaks at 7, 14, 21
#b
inpos.diff = artrans.wge(us$positiveIncrease, phi.tr = 1, lag.max = 100)

c. Seaonality Transformation: stationary data, and an ACF that reflects data other than whitenoise to be modeled.

#c
inpos.diff.seas = artrans.wge(inpos.diff,phi.tr = c(0,0,0,0,0,0,1), lag.max = 100)

d. Model IDing: diagnose with aic.wge to determine phi, thetas: Both AIC/BIC select ARMA(4,2)

#d
aic5.wge(inpos.diff.seas)
## ---------WORKING... PLEASE WAIT... 
## 
## 
## Five Smallest Values of  aic
##       p    q        aic
## 15    4    2   15.62368
## 9     2    2   15.74701
## 18    5    2   15.74771
## 12    3    2   15.75351
## 7     2    0   15.75551
aic5.wge(inpos.diff.seas, type = "bic")
## ---------WORKING... PLEASE WAIT... 
## 
## 
## Five Smallest Values of  bic
##       p    q        bic
## 15    4    2   15.74832
## 2     0    1   15.79201
## 7     2    0   15.80893
## 4     1    0   15.81457
## 3     0    2   15.81522
e. White Noise Test: Reject the H_null and accept this data is not white noise
#e
acf(inpos.diff.seas)

ljung.wge(inpos.diff.seas)$pval
## Obs -0.3711871 -0.02127283 0.07854862 0.05194601 -0.1044838 0.3276797 -0.4658847 0.1565306 0.1284186 -0.1060791 -0.07693327 0.1030723 -0.06391648 0.05002847 0.004471288 -0.03640442 -0.01551996 0.04133725 -0.01110058 -0.08153895 0.06570572 -0.03709739 -0.02390062 0.04588417
## [1] 1.246447e-12
f. Estiamte Phis, Thetas
g. Model Building
#f
est.inpos.seas = est.arma.wge(inpos.diff.seas, p = 4, q = 2)
## 
## Coefficients of Original polynomial:  
## -1.8727 -1.4974 -0.4735 0.0283 
## 
## Factor                 Roots                Abs Recip    System Freq 
## 1+1.1293B+0.6983B^2   -0.8086+-0.8822i      0.8357       0.3681
## 1+0.7944B             -1.2588               0.7944       0.5000
## 1-0.0510B              19.6244               0.0510       0.0000
##   
## 
mean(us$positiveIncrease)
## [1] 22567.12
#g
h. Forecast 7/90 Days 
#7 day
inpos.preds7 = fore.aruma.wge(us$positiveIncrease, phi = est.inpos.seas$phi, theta = est.inpos.seas$theta, d=1, s=7, n.ahead = 7)

#90 day
inpos.preds90 = fore.aruma.wge(us$positiveIncrease, phi = est.inpos.seas$phi, theta = est.inpos.seas$theta, d=1, s=7, n.ahead = 90)

i. Plotting forecasts

#7 day forecast
plot(seq(1,187,1), us$positiveIncrease, type = "l", xlim = c(0,195), ylim = c(0,80000), ylab = "Increase in COVID Cases", main = "7 Day Increase in COVID Cases Forecast")
lines(seq(188, 194,1), inpos.preds7$f, type = "l", col = "red")

#90 day forecast
plot(seq(1,187,1), us$positiveIncrease, type = "l", xlim = c(0,280), ylim = c(0,80000), ylab = "Increase in COVID Cases", main = "90 Day Increase in COVID Cases Forecast")
lines(seq(188, 277,1), inpos.preds90$f, type = "l", col = "red")

Forecast Currently Hospitalized Patients using Positive Increase and Total Cases

  1. Data prep and recognizing differencing and seaonal components of Currently Hospitalized
#Selecting only those dates with reported current hospitilizations
us <- dplyr::slice(us,56:187)
invisible(us)
#Stationarize Currently Hospitalized Variable
us.diff = artrans.wge(us$hospitalizedCurrently, phi.tr = 1, lag.max = 100)

acf(us.diff)
us.diff.seas = artrans.wge(us.diff,phi.tr = c(0,0,0,0,0,0,1), lag.max = 100)

acf(us.diff.seas)
#Stationarize Increase Positive Variable
inpos.diff = artrans.wge(us$positiveIncrease, phi.tr = 1, lag.max = 100)

acf(inpos.diff)
inpos.diff.seas = artrans.wge(inpos.diff,phi.tr = c(0,0,0,0,0,0,1), lag.max = 100)

acf(inpos.diff.seas)

  1. Fit simple regression model predicting hospitalizedCurrently with positiveIncrease.
mlr.fit = lm(us.diff.seas~inpos.diff.seas, data=cbind(data.frame(us.diff.seas), data.frame(inpos.diff.seas)))
summary(mlr.fit)
## 
## Call:
## lm(formula = us.diff.seas ~ inpos.diff.seas, data = cbind(data.frame(us.diff.seas), 
##     data.frame(inpos.diff.seas)))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7149.3  -554.2    73.4   611.2  6506.7 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)   
## (Intercept)     -13.32614  143.23299  -0.093  0.92603   
## inpos.diff.seas   0.11733    0.04227   2.776  0.00637 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1595 on 122 degrees of freedom
## Multiple R-squared:  0.0594, Adjusted R-squared:  0.0517 
## F-statistic: 7.705 on 1 and 122 DF,  p-value: 0.006375
AIC(mlr.fit)
## [1] 2184.712
acf(mlr.fit$residuals)

plot(mlr.fit$residuals)

  1. Diagnose with aic.wge and see what the p, q would be with residuals from above model
    1. Below we can see that our aic.wge function has selected an ARMA(1,2) model for modeling our cmort information.
mlr.phis= aic.wge(mlr.fit$residuals)
mlr.phis
## $type
## [1] "aic"
## 
## $value
## [1] 14.53403
## 
## $p
## [1] 5
## 
## $q
## [1] 2
## 
## $phi
## [1] -0.2465069 -0.4219100  0.3119255  0.2116968  0.2024803
## 
## $theta
## [1] -0.4657537 -0.9566014
## 
## $vara
## [1] 1803071
  1. Now forcast with arima function with phi’s from above coefficients and ARMA(1,2) model,
mlr.fit.arima = arima(us$hospitalizedCurrently, order = c(5,1,2), seasonal = list(order = c(1,0,0), period = 7), xreg = us$positiveIncrease)
## Warning in log(s2): NaNs produced
AIC(mlr.fit.arima)
## [1] 2224.09
  1. Run test to see if data is white noise: confirmed white noise, continue forward
acf(mlr.fit.arima$resid) 

ltest1 = ljung.wge(mlr.fit.arima$resid) 
## Obs -0.006649913 -0.009334226 -0.01070918 0.01827832 -0.02771128 -0.01061308 -0.01123975 0.08781402 -0.05809279 0.07500195 0.05661161 -0.08527475 -0.007553036 0.1098068 0.002427647 -0.010962 -0.04242877 -0.001526285 -0.05612271 -0.03313627 0.02930685 -0.007307643 -0.06056842 -0.03342384
ltest1$pval
## [1] 0.999208
ltest2 = ljung.wge(mlr.fit.arima$resid, K= 48)
## Obs -0.006649913 -0.009334226 -0.01070918 0.01827832 -0.02771128 -0.01061308 -0.01123975 0.08781402 -0.05809279 0.07500195 0.05661161 -0.08527475 -0.007553036 0.1098068 0.002427647 -0.010962 -0.04242877 -0.001526285 -0.05612271 -0.03313627 0.02930685 -0.007307643 -0.06056842 -0.03342384 0.08407797 -0.006875273 -0.002312484 0.06777622 -0.09591675 -0.06198868 -0.01331578 -0.03441461 -0.04091228 0.0270524 0.06615092 0.06538511 -0.09261011 -0.04939741 -0.03187124 -0.04934637 -0.01834029 -0.02031817 -0.001873916 0.02337696 -0.004828854 0.006148926 -0.002455308 0.004810674
ltest2$pval
## [1] 0.9999866
  1. Load the forecasted increase positive in a data frame
    1. 7 Day
#7 Day Case Increase
regs7 = data.frame(positiveIncrease = inpos.preds7$f)
invisible(regs7)
b. 90 Day
#90 Day Case Increase
regs90 = data.frame(positiveIncrease = inpos.preds90$f)
invisible(regs90)
  1. Predictions
    1. 7 Day
mlr1.preds7 = predict(mlr.fit.arima, newxreg = regs7, n.ahead = 7)
invisible(mlr1.preds7)
b. 90 Day
mlr1.preds90 = predict(mlr.fit.arima, newxreg = regs90, n.ahead =90)
invisible(mlr1.preds90)
  1. Plotted Forecasts
    1. 7 Day
plot(seq(1,132,1), us$hospitalizedCurrently, type = "l",  xlim = c(0,140), ylim = c(0,60000), ylab = "Currently Hospitalized COVID Cases", main = "7 Day Forecast for COVID Hospitalized Cases")
lines(seq(133,139,1), mlr1.preds7$pred, type = "l", col = "red")

b. 90 Day

plot(seq(1,132,1), us$hospitalizedCurrently, type = "l",  xlim = c(0,223), ylim = c(0,60000), ylab = "Currently Hospitalized COVID Cases", main = "90 Day Forecast for COVID Hospitalized Cases")
lines(seq(133,222,1), mlr1.preds90$pred, type = "l", col = "red")

  1. Windowed ASE

MLR w/ Correlated Errors for Currently Hospitalized Patients w/ Trend

  1. Fit simple regression model predicting hospitalizedCurrently with positiveIncrease and trend
time <- seq(1,124,1)
mlr.fit.t = lm(us.diff.seas~inpos.diff.seas+time, data=cbind(data.frame(us.diff.seas), data.frame(inpos.diff.seas)))
summary(mlr.fit.t)
## 
## Call:
## lm(formula = us.diff.seas ~ inpos.diff.seas + time, data = cbind(data.frame(us.diff.seas), 
##     data.frame(inpos.diff.seas)))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7111.0  -524.5    73.9   617.7  6540.6 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)   
## (Intercept)      32.35112  289.27171   0.112  0.91114   
## inpos.diff.seas   0.11724    0.04244   2.762  0.00663 **
## time             -0.73096    4.01658  -0.182  0.85590   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1601 on 121 degrees of freedom
## Multiple R-squared:  0.05966,    Adjusted R-squared:  0.04412 
## F-statistic: 3.839 on 2 and 121 DF,  p-value: 0.02419
AIC(mlr.fit.t)
## [1] 2186.678
acf(mlr.fit.t$residuals)

plot(mlr.fit.t$residuals)

  1. Diagnose with aic.wge and see what the p, q would be with residuals from above model
    1. Below we can see that our aic.wge function has selected an ARMA(1,2) model for modeling our cmort information.
mlr.phis= aic.wge(mlr.fit.t$residuals)
mlr.phis
## $type
## [1] "aic"
## 
## $value
## [1] 14.53329
## 
## $p
## [1] 5
## 
## $q
## [1] 2
## 
## $phi
## [1] -0.2466126 -0.4220817  0.3118538  0.2119376  0.2027949
## 
## $theta
## [1] -0.4655288 -0.9564592
## 
## $vara
## [1] 1801724
  1. Now forcast with arima function with phi’s from above coefficients and ARMA(1,2) model,
Time <- seq(1,132,1)
mlr.fit.arima.t = arima(us$hospitalizedCurrently, order = c(5,1,2), seasonal = list(order = c(1,0,0), period = 7), xreg = cbind(us$positiveIncrease, Time))
AIC(mlr.fit.arima.t)
## [1] 2230.756
  1. Run test to see if data is white noise; confirmed white noise continue forward.
acf(mlr.fit.arima.t$resid) 

ltest1 = ljung.wge(mlr.fit.arima.t$resid) 
## Obs 0.002267493 0.004743329 0.004085712 0.008584111 0.003613306 0.03592544 -0.02552583 0.08439651 -0.1060931 0.01276058 0.02640347 -0.06061388 0.05268146 0.1458814 0.008052735 -0.04255608 -0.0935536 -0.04208266 -0.0498013 -0.001717023 0.05934226 0.01341049 -0.07137289 -0.07279666
ltest1$pval
## [1] 0.9823161
ltest2 = ljung.wge(mlr.fit.arima.t$resid, K= 48)
## Obs 0.002267493 0.004743329 0.004085712 0.008584111 0.003613306 0.03592544 -0.02552583 0.08439651 -0.1060931 0.01276058 0.02640347 -0.06061388 0.05268146 0.1458814 0.008052735 -0.04255608 -0.0935536 -0.04208266 -0.0498013 -0.001717023 0.05934226 0.01341049 -0.07137289 -0.07279666 0.05320601 -0.01313386 0.01794322 0.1027134 -0.07595104 -0.06199621 -0.03306343 -0.06206906 -0.03944973 0.03988941 0.07572305 0.07297186 -0.1029145 -0.06857141 -0.04403521 -0.04289784 0.005428483 -0.003871671 -7.833593e-05 0.005208028 -0.03085214 -0.01036108 0.00645461 0.02868683
ltest2$pval
## [1] 0.9990337
  1. Load the forecasted increase positive in a data frame
    1. 7 Day
#7 Day Case Increase
regs7t = data.frame(cbind(positiveIncrease = inpos.preds7$f, Time = seq(133,139)))
invisible(regs7)
b. 90 Day
#90 Day Case Increase
regs90t = data.frame(cbind(positiveIncrease = inpos.preds90$f, Time = seq(133,222)))
invisible(regs90)
  1. Predictions
    1. 7 Day
mlr1.preds7.t = predict(mlr.fit.arima.t, newxreg = regs7t, n.ahead = 7)
invisible(mlr1.preds7.t)
b. 90 Day
mlr1.preds90.t = predict(mlr.fit.arima.t, newxreg = regs90t, n.ahead =90)
invisible(mlr1.preds90.t)
  1. Plotted Forecasts
    1. 7 Day
plot(seq(1,132,1), us$hospitalizedCurrently, type = "l",  xlim = c(0,140), ylim = c(0,60000), ylab = "Currently Hospitalized COVID Cases", main = "7 Day Forecast for COVID Hospitalized Cases")
lines(seq(133,139,1), mlr1.preds7.t$pred, type = "l", col = "red")

b. 90 Day

plot(seq(1,132,1), us$hospitalizedCurrently, type = "l",  xlim = c(0,223), ylim = c(0,85000), ylab = "Currently Hospitalized COVID Cases", main = "90 Day Forecast for COVID Hospitalized Cases")
lines(seq(133,222,1), mlr1.preds90.t$pred, type = "l", col = "red")

MLR w/ Correlated Errors (lagged variables) for Currently Hospitalized Patients

With a quick check, we can see that there is no lag correlationed between the increase of COVID patients and hospitalized patients, like we theorized there might be. So we will not model an MLR w/ correlated errors on lagged variables.

ccf(us$positiveIncrease,us$hospitalizedCurrently)

VAR Model

Since we are working with a seasonal and transformation component but it requires an additional transformation for the total positive COVID cases to become stationary for evaluation, we will only use positive increase of cases to predict currently hospitzliaed cases for the VAR model.

  1. Differenced and Seasonal Transformation VAR Model
#Positive Cases Transformations
#pos.diff = artrans.wge(us$positive, phi.tr = 1, lag.max = 100)
#pos.diff.2 = artrans.wge(pos.diff, phi.tr = 1, lag.max = 100)
#pos.trans = artrans.wge(pos.diff.2,phi.tr = c(0,0,0,0,0,0,1), lag.max = 100)
#Increase Positive Transformation
inpos.diff = artrans.wge(us$positiveIncrease, phi.tr = 1, lag.max = 100)

inpos.trans = artrans.wge(inpos.diff,phi.tr = c(0,0,0,0,0,0,1), lag.max = 100)

#Current Hospitalization Transformations
us.diff = artrans.wge(us$hospitalizedCurrently, phi.tr = 1, lag.max = 100)

currhosp.trans = artrans.wge(us.diff,phi.tr = c(0,0,0,0,0,0,1), lag.max = 100)

  1. Use Var select to find best model and fit model: p = 8 for lowest AIC
#VARselect to choose lag
VARselect(cbind(currhosp.trans, inpos.trans), lag.max = 10, type= "both")
## $selection
## AIC(n)  HQ(n)  SC(n) FPE(n) 
##      7      7      2      7 
## 
## $criteria
##                   1            2            3            4            5
## AIC(n) 3.093857e+01 3.083329e+01 3.079091e+01 3.076976e+01 3.069554e+01
## HQ(n)  3.101650e+01 3.095018e+01 3.094676e+01 3.096458e+01 3.092932e+01
## SC(n)  3.113058e+01 3.112131e+01 3.117493e+01 3.124979e+01 3.127158e+01
## FPE(n) 2.731960e+13 2.459304e+13 2.357880e+13 2.309552e+13 2.145770e+13
##                   6            7            8            9           10
## AIC(n) 3.066656e+01 3.051427e+01 3.053424e+01 3.059674e+01 3.059649e+01
## HQ(n)  3.093930e+01 3.082598e+01 3.088492e+01 3.098638e+01 3.102509e+01
## SC(n)  3.133860e+01 3.128233e+01 3.139830e+01 3.155681e+01 3.165257e+01
## FPE(n) 2.086403e+13 1.793906e+13 1.833019e+13 1.955149e+13 1.959499e+13
#fit model
usdiff.var.fit = VAR(cbind(currhosp.trans, inpos.trans), type = "both", p = 2)
  1. Predictions for Difference
    1. 7 Day
#7 day predictions
pred.var7 = predict(usdiff.var.fit, n.ahead = 7)
pred.var7$fcst$currhosp.trans
##            fcst     lower    upper       CI
## [1,]  -42.12350 -2942.410 2858.163 2900.286
## [2,] -385.17371 -3346.493 2576.145 2961.319
## [3,]  -80.65766 -3262.529 3101.214 3181.871
## [4,] -124.14161 -3327.324 3079.041 3203.182
## [5,]  -46.62794 -3284.538 3191.282 3237.910
## [6,]  -43.58380 -3288.207 3201.039 3244.623
## [7,]  -11.23201 -3262.277 3239.813 3251.045
b. 90 Day
#90 day predictions
pred.var90 = predict(usdiff.var.fit, n.ahead = 90)
pred.var90$fcst$currhosp.trans
##              fcst     lower    upper       CI
##  [1,]  -42.123501 -2942.410 2858.163 2900.286
##  [2,] -385.173707 -3346.493 2576.145 2961.319
##  [3,]  -80.657665 -3262.529 3101.214 3181.871
##  [4,] -124.141611 -3327.324 3079.041 3203.182
##  [5,]  -46.627942 -3284.538 3191.282 3237.910
##  [6,]  -43.583795 -3288.207 3201.039 3244.623
##  [7,]  -11.232008 -3262.277 3239.813 3251.045
##  [8,]   -6.349200 -3259.226 3246.528 3252.877
##  [9,]    7.305513 -3246.870 3261.481 3254.175
## [10,]   12.421029 -3242.210 3267.053 3254.632
## [11,]   18.773716 -3236.131 3273.679 3254.905
## [12,]   22.480950 -3232.533 3277.495 3255.014
## [13,]   26.203504 -3228.870 3281.277 3255.073
## [14,]   28.923414 -3226.175 3284.022 3255.099
## [15,]   31.503265 -3223.608 3286.615 3255.112
## [16,]   33.699320 -3221.418 3288.817 3255.117
## [17,]   35.776726 -3219.344 3290.897 3255.120
## [18,]   37.697587 -3217.424 3292.819 3255.122
## [19,]   39.549962 -3215.572 3294.672 3255.122
## [20,]   41.334045 -3213.789 3296.457 3255.123
## [21,]   43.082251 -3212.041 3298.205 3255.123
## [22,]   44.799892 -3210.323 3299.923 3255.123
## [23,]   46.499486 -3208.623 3301.622 3255.123
## [24,]   48.185069 -3206.938 3303.308 3255.123
## [25,]   49.861824 -3205.261 3304.985 3255.123
## [26,]   51.532055 -3203.591 3306.655 3255.123
## [27,]   53.198022 -3201.925 3308.321 3255.123
## [28,]   54.860928 -3200.262 3309.984 3255.123
## [29,]   56.521787 -3198.601 3311.645 3255.123
## [30,]   58.181203 -3196.942 3313.304 3255.123
## [31,]   59.839642 -3195.283 3314.963 3255.123
## [32,]   61.497398 -3193.626 3316.620 3255.123
## [33,]   63.154688 -3191.968 3318.278 3255.123
## [34,]   64.811654 -3190.311 3319.935 3255.123
## [35,]   66.468399 -3188.655 3321.591 3255.123
## [36,]   68.124990 -3186.998 3323.248 3255.123
## [37,]   69.781476 -3185.341 3324.904 3255.123
## [38,]   71.437889 -3183.685 3326.561 3255.123
## [39,]   73.094252 -3182.029 3328.217 3255.123
## [40,]   74.750580 -3180.372 3329.874 3255.123
## [41,]   76.406884 -3178.716 3331.530 3255.123
## [42,]   78.063172 -3177.060 3333.186 3255.123
## [43,]   79.719449 -3175.404 3334.842 3255.123
## [44,]   81.375717 -3173.747 3336.499 3255.123
## [45,]   83.031981 -3172.091 3338.155 3255.123
## [46,]   84.688240 -3170.435 3339.811 3255.123
## [47,]   86.344498 -3168.778 3341.467 3255.123
## [48,]   88.000753 -3167.122 3343.124 3255.123
## [49,]   89.657007 -3165.466 3344.780 3255.123
## [50,]   91.313260 -3163.810 3346.436 3255.123
## [51,]   92.969513 -3162.153 3348.092 3255.123
## [52,]   94.625766 -3160.497 3349.749 3255.123
## [53,]   96.282018 -3158.841 3351.405 3255.123
## [54,]   97.938270 -3157.185 3353.061 3255.123
## [55,]   99.594521 -3155.528 3354.717 3255.123
## [56,]  101.250773 -3153.872 3356.374 3255.123
## [57,]  102.907025 -3152.216 3358.030 3255.123
## [58,]  104.563276 -3150.560 3359.686 3255.123
## [59,]  106.219528 -3148.903 3361.342 3255.123
## [60,]  107.875779 -3147.247 3362.999 3255.123
## [61,]  109.532031 -3145.591 3364.655 3255.123
## [62,]  111.188282 -3143.935 3366.311 3255.123
## [63,]  112.844534 -3142.278 3367.967 3255.123
## [64,]  114.500785 -3140.622 3369.624 3255.123
## [65,]  116.157037 -3138.966 3371.280 3255.123
## [66,]  117.813288 -3137.310 3372.936 3255.123
## [67,]  119.469540 -3135.653 3374.592 3255.123
## [68,]  121.125791 -3133.997 3376.249 3255.123
## [69,]  122.782043 -3132.341 3377.905 3255.123
## [70,]  124.438294 -3130.685 3379.561 3255.123
## [71,]  126.094546 -3129.028 3381.218 3255.123
## [72,]  127.750797 -3127.372 3382.874 3255.123
## [73,]  129.407049 -3125.716 3384.530 3255.123
## [74,]  131.063300 -3124.060 3386.186 3255.123
## [75,]  132.719552 -3122.403 3387.843 3255.123
## [76,]  134.375803 -3120.747 3389.499 3255.123
## [77,]  136.032055 -3119.091 3391.155 3255.123
## [78,]  137.688306 -3117.435 3392.811 3255.123
## [79,]  139.344558 -3115.778 3394.468 3255.123
## [80,]  141.000809 -3114.122 3396.124 3255.123
## [81,]  142.657061 -3112.466 3397.780 3255.123
## [82,]  144.313312 -3110.810 3399.436 3255.123
## [83,]  145.969564 -3109.153 3401.093 3255.123
## [84,]  147.625815 -3107.497 3402.749 3255.123
## [85,]  149.282067 -3105.841 3404.405 3255.123
## [86,]  150.938318 -3104.185 3406.061 3255.123
## [87,]  152.594570 -3102.528 3407.718 3255.123
## [88,]  154.250821 -3100.872 3409.374 3255.123
## [89,]  155.907073 -3099.216 3411.030 3255.123
## [90,]  157.563324 -3097.560 3412.686 3255.123
  1. Calculate Actual Forecasts from Predicted Differences
    1. 7 Day
startingPoints7 = us$hospitalizedCurrently[126:132]
currHospForecasts7 = pred.var7$fcst$currhosp.trans[,1:3] + startingPoints7
b. 90 Day
startingPoints90 = us$hospitalizedCurrently[43:132]
currHospForecasts90 = pred.var90$fcst$currhosp.trans[,1:3] + startingPoints90
  1. Plotting Forecasts
    1. 7 Day
#7 day Forecasts
plot(seq(1,132,1), us$hospitalizedCurrently, type = "l", xlim = c(0,139), ylab = "Currently Hospitalized COVID Patients", main = "7 Day Currently Hospitalized Patients Forecast")
lines(seq(133,139,1), as.data.frame(currHospForecasts7)$fcst, type = "l", col = "red")

b. 90 Day

#90 day Forecasts
plot(seq(1,132,1), us$hospitalizedCurrently, type = "l", xlim = c(0,222), ylab = "Currently Hospitalized COVID Patients", main = "90 Day Currently Hospitalized Patients Forecast")
lines(seq(133,222,1), as.data.frame(currHospForecasts90)$fcst, type = "l", col = "red")

  1. Calculate ASE # WHERE YOU LEFT OFF
#varASE = mean((us$hospitalizedCurrently[123:132]-pred.var$fcst$y1[1:10])^2)
#varASE
  1. Windowed ASE

Multilayer Perceptron Model (MLP) / NN

Since we have been looking at additional regressors for all our above models and know that the best regressor to leverage will be positive increase of COVID cases. This regressor reflects similar behavior to that of current hospitlizations and can be model properly against hospitalizations. 1. Create Current Hospitalized ts variable

tUScurrhop = ts(us$hospitalizedCurrently)
  1. Create data frame of regressor: positive increase covid cases variable
tUSx = data.frame(positiveIncrease = ts(us$positiveIncrease))
  1. Forecast of positive increase of COVID cases
    1. 7 Day Forecast for Regressor
fit.mlp.new = mlp(ts(tUSx), reps = 50, comb = "mean")
mlp.new.fore7 = forecast(fit.mlp.new, h = 7)
mlp.new.fore7
##     Point Forecast
## 133       60120.38
## 134       63466.16
## 135       65066.40
## 136       66008.94
## 137       65216.75
## 138       64832.34
## 139       64545.57
b. 90 Day Forecast for Regressor
mlp.new.fore90 = forecast(fit.mlp.new, h = 90)
mlp.new.fore90
##     Point Forecast
## 133       60120.38
## 134       63466.16
## 135       65066.40
## 136       66008.94
## 137       65216.75
## 138       64832.34
## 139       64545.57
## 140       64840.99
## 141       64982.79
## 142       65111.03
## 143       65029.59
## 144       65001.26
## 145       64965.07
## 146       65001.73
## 147       65007.93
## 148       65018.30
## 149       65000.98
## 150       65000.58
## 151       64998.12
## 152       65006.27
## 153       65004.16
## 154       65003.08
## 155       64997.89
## 156       64999.09
## 157       64999.73
## 158       65001.94
## 159       65000.23
## 160       64998.96
## 161       64997.27
## 162       64998.12
## 163       64998.69
## 164       64999.35
## 165       64998.45
## 166       64997.80
## 167       64997.23
## 168       64997.66
## 169       64997.93
## 170       64998.12
## 171       64997.71
## 172       64997.45
## 173       64997.26
## 174       64997.46
## 175       64997.56
## 176       64997.60
## 177       64997.43
## 178       64997.33
## 179       64997.28
## 180       64997.37
## 181       64997.39
## 182       64997.40
## 183       64997.32
## 184       64997.29
## 185       64997.28
## 186       64997.32
## 187       64997.32
## 188       64997.32
## 189       64997.29
## 190       64997.28
## 191       64997.28
## 192       64997.29
## 193       64997.30
## 194       64997.29
## 195       64997.28
## 196       64997.28
## 197       64997.28
## 198       64997.28
## 199       64997.28
## 200       64997.28
## 201       64997.28
## 202       64997.28
## 203       64997.28
## 204       64997.28
## 205       64997.28
## 206       64997.28
## 207       64997.27
## 208       64997.28
## 209       64997.28
## 210       64997.28
## 211       64997.28
## 212       64997.28
## 213       64997.27
## 214       64997.27
## 215       64997.28
## 216       64997.28
## 217       64997.28
## 218       64997.27
## 219       64997.27
## 220       64997.27
## 221       64997.27
## 222       64997.27
  1. Combine observed new cases + forecast new cases
    1. 7 Day regressor var
new.regressor7 <- data.frame(c(us$positiveIncrease, mlp.new.fore7$mean))
new.regressor7 
##     c.us.positiveIncrease..mlp.new.fore7.mean.
## 1                                      3613.00
## 2                                      3171.00
## 3                                      4671.00
## 4                                      6255.00
## 5                                      6885.00
## 6                                      9259.00
## 7                                     11442.00
## 8                                     10632.00
## 9                                     12873.00
## 10                                    17656.00
## 11                                    19044.00
## 12                                    19724.00
## 13                                    19655.00
## 14                                    21897.00
## 15                                    24694.00
## 16                                    25773.00
## 17                                    27992.00
## 18                                    31945.00
## 19                                    33252.00
## 20                                    25550.00
## 21                                    28937.00
## 22                                    30737.00
## 23                                    30534.00
## 24                                    34419.00
## 25                                    34351.00
## 26                                    30561.00
## 27                                    27844.00
## 28                                    25133.00
## 29                                    25631.00
## 30                                    30298.00
## 31                                    30923.00
## 32                                    31964.00
## 33                                    27963.00
## 34                                    27468.00
## 35                                    25872.00
## 36                                    26333.00
## 37                                    28916.00
## 38                                    31784.00
## 39                                    34174.00
## 40                                    35901.00
## 41                                    27380.00
## 42                                    22605.00
## 43                                    25219.00
## 44                                    26641.00
## 45                                    29568.00
## 46                                    33056.00
## 47                                    29185.00
## 48                                    25767.00
## 49                                    22523.00
## 50                                    22449.00
## 51                                    25056.00
## 52                                    27490.00
## 53                                    27605.00
## 54                                    24810.00
## 55                                    21504.00
## 56                                    18236.00
## 57                                    22663.00
## 58                                    21193.00
## 59                                    26705.00
## 60                                    24710.00
## 61                                    24701.00
## 62                                    20171.00
## 63                                    20992.00
## 64                                    20853.00
## 65                                    21319.00
## 66                                    26513.00
## 67                                    24555.00
## 68                                    21590.00
## 69                                    20105.00
## 70                                    18798.00
## 71                                    16629.00
## 72                                    19385.00
## 73                                    22601.00
## 74                                    23522.00
## 75                                    23762.00
## 76                                    21639.00
## 77                                    20415.00
## 78                                    19982.00
## 79                                    20312.00
## 80                                    20839.00
## 81                                    23352.00
## 82                                    23106.00
## 83                                    18823.00
## 84                                    17012.00
## 85                                    17175.00
## 86                                    20803.00
## 87                                    22032.00
## 88                                    23477.00
## 89                                    25341.00
## 90                                    21373.00
## 91                                    18510.00
## 92                                    23435.00
## 93                                    23873.00
## 94                                    27527.00
## 95                                    31046.00
## 96                                    31994.00
## 97                                    27284.00
## 98                                    27017.00
## 99                                    33021.00
## 100                                   38684.00
## 101                                   39072.00
## 102                                   44421.00
## 103                                   43783.00
## 104                                   41857.00
## 105                                   39175.00
## 106                                   47462.00
## 107                                   50674.00
## 108                                   53826.00
## 109                                   54223.00
## 110                                   54734.00
## 111                                   45789.00
## 112                                   41600.00
## 113                                   51766.00
## 114                                   62147.00
## 115                                   58836.00
## 116                                   66645.00
## 117                                   63007.00
## 118                                   60978.00
## 119                                   58465.00
## 120                                   62879.00
## 121                                   65382.00
## 122                                   70953.00
## 123                                   77233.00
## 124                                   65180.00
## 125                                   64884.00
## 126                                   56971.00
## 127                                   63642.00
## 128                                   69150.00
## 129                                   71027.00
## 130                                   75193.00
## 131                                   65413.00
## 132                                   61713.00
## 133                                   60120.38
## 134                                   63466.16
## 135                                   65066.40
## 136                                   66008.94
## 137                                   65216.75
## 138                                   64832.34
## 139                                   64545.57
b. 90 day regressor var
new.regressor90 <- data.frame(c(us$positiveIncrease, mlp.new.fore90$mean))
new.regressor90
##     c.us.positiveIncrease..mlp.new.fore90.mean.
## 1                                       3613.00
## 2                                       3171.00
## 3                                       4671.00
## 4                                       6255.00
## 5                                       6885.00
## 6                                       9259.00
## 7                                      11442.00
## 8                                      10632.00
## 9                                      12873.00
## 10                                     17656.00
## 11                                     19044.00
## 12                                     19724.00
## 13                                     19655.00
## 14                                     21897.00
## 15                                     24694.00
## 16                                     25773.00
## 17                                     27992.00
## 18                                     31945.00
## 19                                     33252.00
## 20                                     25550.00
## 21                                     28937.00
## 22                                     30737.00
## 23                                     30534.00
## 24                                     34419.00
## 25                                     34351.00
## 26                                     30561.00
## 27                                     27844.00
## 28                                     25133.00
## 29                                     25631.00
## 30                                     30298.00
## 31                                     30923.00
## 32                                     31964.00
## 33                                     27963.00
## 34                                     27468.00
## 35                                     25872.00
## 36                                     26333.00
## 37                                     28916.00
## 38                                     31784.00
## 39                                     34174.00
## 40                                     35901.00
## 41                                     27380.00
## 42                                     22605.00
## 43                                     25219.00
## 44                                     26641.00
## 45                                     29568.00
## 46                                     33056.00
## 47                                     29185.00
## 48                                     25767.00
## 49                                     22523.00
## 50                                     22449.00
## 51                                     25056.00
## 52                                     27490.00
## 53                                     27605.00
## 54                                     24810.00
## 55                                     21504.00
## 56                                     18236.00
## 57                                     22663.00
## 58                                     21193.00
## 59                                     26705.00
## 60                                     24710.00
## 61                                     24701.00
## 62                                     20171.00
## 63                                     20992.00
## 64                                     20853.00
## 65                                     21319.00
## 66                                     26513.00
## 67                                     24555.00
## 68                                     21590.00
## 69                                     20105.00
## 70                                     18798.00
## 71                                     16629.00
## 72                                     19385.00
## 73                                     22601.00
## 74                                     23522.00
## 75                                     23762.00
## 76                                     21639.00
## 77                                     20415.00
## 78                                     19982.00
## 79                                     20312.00
## 80                                     20839.00
## 81                                     23352.00
## 82                                     23106.00
## 83                                     18823.00
## 84                                     17012.00
## 85                                     17175.00
## 86                                     20803.00
## 87                                     22032.00
## 88                                     23477.00
## 89                                     25341.00
## 90                                     21373.00
## 91                                     18510.00
## 92                                     23435.00
## 93                                     23873.00
## 94                                     27527.00
## 95                                     31046.00
## 96                                     31994.00
## 97                                     27284.00
## 98                                     27017.00
## 99                                     33021.00
## 100                                    38684.00
## 101                                    39072.00
## 102                                    44421.00
## 103                                    43783.00
## 104                                    41857.00
## 105                                    39175.00
## 106                                    47462.00
## 107                                    50674.00
## 108                                    53826.00
## 109                                    54223.00
## 110                                    54734.00
## 111                                    45789.00
## 112                                    41600.00
## 113                                    51766.00
## 114                                    62147.00
## 115                                    58836.00
## 116                                    66645.00
## 117                                    63007.00
## 118                                    60978.00
## 119                                    58465.00
## 120                                    62879.00
## 121                                    65382.00
## 122                                    70953.00
## 123                                    77233.00
## 124                                    65180.00
## 125                                    64884.00
## 126                                    56971.00
## 127                                    63642.00
## 128                                    69150.00
## 129                                    71027.00
## 130                                    75193.00
## 131                                    65413.00
## 132                                    61713.00
## 133                                    60120.38
## 134                                    63466.16
## 135                                    65066.40
## 136                                    66008.94
## 137                                    65216.75
## 138                                    64832.34
## 139                                    64545.57
## 140                                    64840.99
## 141                                    64982.79
## 142                                    65111.03
## 143                                    65029.59
## 144                                    65001.26
## 145                                    64965.07
## 146                                    65001.73
## 147                                    65007.93
## 148                                    65018.30
## 149                                    65000.98
## 150                                    65000.58
## 151                                    64998.12
## 152                                    65006.27
## 153                                    65004.16
## 154                                    65003.08
## 155                                    64997.89
## 156                                    64999.09
## 157                                    64999.73
## 158                                    65001.94
## 159                                    65000.23
## 160                                    64998.96
## 161                                    64997.27
## 162                                    64998.12
## 163                                    64998.69
## 164                                    64999.35
## 165                                    64998.45
## 166                                    64997.80
## 167                                    64997.23
## 168                                    64997.66
## 169                                    64997.93
## 170                                    64998.12
## 171                                    64997.71
## 172                                    64997.45
## 173                                    64997.26
## 174                                    64997.46
## 175                                    64997.56
## 176                                    64997.60
## 177                                    64997.43
## 178                                    64997.33
## 179                                    64997.28
## 180                                    64997.37
## 181                                    64997.39
## 182                                    64997.40
## 183                                    64997.32
## 184                                    64997.29
## 185                                    64997.28
## 186                                    64997.32
## 187                                    64997.32
## 188                                    64997.32
## 189                                    64997.29
## 190                                    64997.28
## 191                                    64997.28
## 192                                    64997.29
## 193                                    64997.30
## 194                                    64997.29
## 195                                    64997.28
## 196                                    64997.28
## 197                                    64997.28
## 198                                    64997.28
## 199                                    64997.28
## 200                                    64997.28
## 201                                    64997.28
## 202                                    64997.28
## 203                                    64997.28
## 204                                    64997.28
## 205                                    64997.28
## 206                                    64997.28
## 207                                    64997.27
## 208                                    64997.28
## 209                                    64997.28
## 210                                    64997.28
## 211                                    64997.28
## 212                                    64997.28
## 213                                    64997.27
## 214                                    64997.27
## 215                                    64997.28
## 216                                    64997.28
## 217                                    64997.28
## 218                                    64997.27
## 219                                    64997.27
## 220                                    64997.27
## 221                                    64997.27
## 222                                    64997.27
  1. Fit model for currently hospitalized w/ regressor of new cases
mlp.fit1 = mlp(tUScurrhop, xreg = tUSx, outplot = T, comb = "mean")

plot(mlp.fit1)

  1. Forecast w/ known Positive Increase Case Variable
    1. Currently Hospitalized 7 Day Forecast w/ Positive Increaser Regressor
currhosp.fore7 = forecast(mlp.fit1, h = 7, xreg = new.regressor7)
plot(currhosp.fore7)

b. Currently Hospitalized 90 Day Forecast w/ Positive Increaser Regressor

currhosp.fore90 = forecast(mlp.fit1, h = 90, xreg = new.regressor90)
plot(currhosp.fore90)

  1. Windowed ASE
  1. 7 Day Forecast
#if (!require(devtools)) {
#  install.packages("devtools")
#}
#devtools::install_github("josephsdavid/tswgewrapped", build_vignettes = TRUE)
library(tswgewrapped)
- Train/ Test Data Sets
data_train.m <- data.frame(hospitalizedCurrently = us$hospitalizedCurrently[1:122], positiveIncrease = us$positiveIncrease[1:122])
data_test.m <- data.frame(hospitalizedCurrently = us$hospitalizedCurrently[123:132], positiveIncrease = us$positiveIncrease[123:132])
# search for best NN hyperparameters in given grid
model.m = tswgewrapped::ModelBuildNNforCaret$new(data = data_train.m, var_interest = "hospitalizedCurrently",
                                               search = 'random', tuneLength = 5, parallel = TRUE,
                                               batch_size = 50, h = 7, m = 7,
                                               verbose = 1)
## Registered S3 method overwritten by 'lava':
##   method     from   
##   print.pcor greybox
## Aggregating results
## Selecting tuning parameters
## Fitting reps = 10, hd = 2, allow.det.season = TRUE on full training set
## - Total Time for training: : 65.979 sec elapsed
- The ASEs associated with the grid of hyperparameters is shown in the table and heatmap below.
res.m <- model.m$summarize_hyperparam_results()
res.m
##   reps hd allow.det.season     RMSE      ASE   RMSESD    ASESD
## 1   10  2             TRUE 3103.429 13970241 2184.689 16296534
## 2   11  3            FALSE 3129.341 14942705 2380.109 19673017
## 3   16  1            FALSE 3559.638 16480775 2047.127 15988977
## 4   16  3             TRUE 3204.752 15391189 2373.358 19509361
## 5   22  3            FALSE 3268.191 16196849 2463.200 19979406
model.m$plot_hyperparam_results()

- Windowed ASE: The best hyperparemeters are shown below
best.m <- model.m$summarize_best_hyperparams()
best.m
##   reps hd allow.det.season
## 1   10  2             TRUE
final.ase.m <- dplyr::filter(res.m, reps == best.m$reps &
                    hd == best.m$hd &
                    allow.det.season == best.m$allow.det.season)[['ASE']]
final.ase.m
## [1] 13970241
- The best hyperparameters based on this grid search are 13 repetitions and 1 hidden layers, and allow.det.season = FALSE which has a mean rolling window ASE of 5.4964555.
# Final Model
caret_model.m = model.m$get_final_models(subset = 'a')
caret_model.m$finalModel
## MLP fit with 2 hidden nodes and 10 repetitions.
## Series modelled in differences: D1.
## Univariate lags: (1,2,3)
## 1 regressor included.
## - Regressor 1 lags: (1,3)
## Forecast combined using the median operator.
## MSE: 1065689.5459.
- Plot final model
# Plot Final Model
plot(caret_model.m$finalModel)